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Abstract 

The three-dimensional four-fermion model with a Z2 chiral symmetry and = 4 fermion 
species is investigated numerically at non-zero temparature T and chemical potential /x. The 
phase diagram in the (/x, T) plane is mapped out quantitatively. A detailed finite size scaling 
analysis shows that the T = 0, 7^ chiral phase transition is first order, whereas the 
T 7^ 0, /i 7^ transition remains second order down to very low T. Quantitative results for 
the location of the end of the first order line are given. The simulations at T = provided 
no evidence for a nuclear liquid-gas transition at < /i^. 
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1 Introduction 



The chiral phase transition in QCD separating high temperature (high density) quark-gluon 
plasma phase from low temperature (low density) hadronic phase has been studied intensively 
in the last decade. At temperatures of the order of the pion mass, T ^ 140 MeV or densities 
a few times that of ordinary nuclear matter, n ~ (3 — 5)nQ with no ~ 0.15 fm~^, conditions 
are reached where hadrons start overlapping, only the short distance interaction among the 
internal parton degrees of freedom described by QCD is important and dynamical breaking 
of chiral symmetry no longer occurs. Understanding the properties of this transition is 
becoming increasingly important in view of recent experimental efforts to create and detect 
the quark-gluon plasma in the relativistic heavy-ion collisions at BNL and CERN. 

Since the problem of chiral symmetry breaking and its restoration is intrinsically non- 
perturbative, the number of available techniques is limited and most of our knowledge about 
the phenomenon comes from lattice simulations. Due to the complexity of QCD, studies have 
so far been done on lattices of modest size and have been unable to yield quantitative results 
as far as the universality class of the finite temperature phase transition is concerned. Also, 
finite density QCD, which is more relevant to real physics, is a field where lattice technology 
is still at an early stage IQ . The complex nature of the determinant of the Dirac operator 
at finite chemical potential makes it impossible to use standard simulation algorithms based 
on positive definite probability functions. 

This manuscript addresses the problem of chiral symmetry restoration at non-zero tem- 
perature and non-zero chemical potential in the three-dimensional four-fermion model 0, Q 
with a Z2 chiral symmetry in order to understand what ingredients might play a decisive 
role in more complex systems like gauge theories. The model has been simplified as much as 
possible in order to produce data of the highest quality and get some insight into the range 
of parameters we need for studies of more realistic problems. The large N description (which 
is equivalent to a mean field approximation) of the phase structure of the three-dimensional 
four-fermion model at non-zero chemical potential predicts a first order transition for 
T = and a continuous transition for T > 0. The critical value of the chemical potential fic 
is equal to the value of the fermion mass at = 0. Interactions as expected decrease /ic below 
the mean field result A recent numerical study of the T 7^ 0, yU = transition favors 
the two-dimensional Ising model exponents for < 00 Q. In this manuscript we present 
the results of Monte Carlo simulations with a small number of fermion species, i.e., A^ = 4 
in order to test the mean field theory predictions and look for possible next-to-the-leading 
order corrections. We show in this paper by employing finite size scaling techniques that the 
T = 0, /i 7^ phase transition is first order. We also show that the second order nature of 
the T 7^ 0, /i = transition remains stable down to low temparatures and large values of 
the chemical potential. A simple continuity argument suggests that a tricritical point in the 
{n, T) plane must lie on the low temperature and high density section of the critical line. 
More specifically, our simulations show that the tricritical point lies on the section of the 
phase boundary defined by: T/Tc < 0.23, fi/fic > 0.97, where Tc is the zero density critical 
temperature and fic is the zero temperature critical chemical potential. Our simulations 
did not provide any evidence for non-zero fermion number density at T = and fi < fic- 
This implies that the nuclear liquid-gas transition is either extremely weak or very close to 



2 



the chiral transition or it does not exist in this model. Our results show that higher order 
corrections in the expansion for the nature and the location of the transition points in 
the phase diagram are very small for this model. 

2 The model 

The continuum spacetime (we work in Euclidean space throughout) Lagrangian of the model 
is 

C = + ll-iQ + mq)i)i - . (1) 

Here ipi is a four- component spinor and the index i runs over N fermion species, /i is the 
chemical potential and mq is the fermion bare mass. For zero fermion bare mass, £, has a 
Z2 chiral symmetry 

-0 ^ 75^ ; ^ ^ -V'Ts (2) 

which is spontaneously broken whenever a non- vanishing condensate {ipip) is generated. Both 
analytical and numerical work in this model are aided by the introduction of an auxiliary 
field 0", so Eq. (|I]) becomes 

£ = ^,(^ + cr)^, + ^a^ (3) 

At tree level, the field a has no dynamics; it is truly an auxiliary field. However, it acquires 
dynamical content by dint of quantum effects arising from integrating out the fermions. In 
the chiral limit the vacuum expectation value of a becomes a dynamical fermion mass and 
it is a convenient order parameter for any chiral phase transition in the theory. 

There are several motivations for studying such a simple model: (i) for sufficiently strong 
coupling it exhibits chiral symmetry breaking at zero temperature and density; (ii) the 
spectrum of excitations contains both baryons and mesons, i.e. the elementary fermions and 
the composite fermion— anti-fermion states; (iii) the model has an interacting continuum 
limit; (iv) when formulated on the lattice it has a real Euclidean action even for non-zero 
chemical potential; (v) the Yukawa-like coupling makes the inverse Dirac operator non- 
singular and the simulations can be done in the chiral limit directly, whereas in gauge theories 
the matrix to be inverted becomes ill-conditioned and the simulations have to be done using a 
fermion bare mass; (vi) the model has fewer degrees of freedom than QCD, and can be studied 
with greater precision on much bigger lattices than the lattices presently used for QCD 
thermodynamics. The chiral symmetry must be discrete in order to observe a phase transition 
at non-zero temperature in ci = 2 + 1, otherwise according to Coleman— Mermin— Wagner 
theorem [|^ a continuous chiral symmetry must be manifested algebraically for all T > 0. 

In four dimensions the four-fermion model is a trivial theory |^ and is believed to be an 
effective theory of quarks and gluons at intermediate energies. Unfortunately, the order of 
the chiral phase transition at non- zero chemical potential in this model depends on the value 
of the cutoff 1^. In two dimensions the discrete chiral symmetry restoration is dominated 
by the materialization of kink— anti-kink states. The zero temperature and finite chemical 
potential transition is first order, whereas for any finite number of fermion species the 
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chiral symmetry is restored at any non-zero temperature due to the condensation of the 



kinks [|10 



The four-fermion model in its bosonized form may be formulated on a spacetime lattice 
using the following action: 



N/2 
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Slat = E E H^)M.,yXi{y) + g E Xi{x)Xi{x) E ^(^) + ^ E ^'(^)' (4) 
where Xi and Xi are Grassmann-valued staggered fermion fields defined on the lattice sites, 



the auxiliary field a is defined on the dual lattice sites |jTT[, and the symbol denotes 
the set of 8 dual lattice sites x surrounding the direct lattice site x. The lattice spacing a 
has been set to one for convenience. Details and motivation of this particular scheme can be 
found in [Q. The fermion kinetic operator M is given by 
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u=1.2 



+ ^ E Vuix) [5y^^ + - 5y^^-0\ , (5) 



where riy{x) are the Kawamoto-Smit phases (^_i^a^o+...+a;^_i_ influence of the chemical 
potential is manifested through the timelinks following |]T2|. Only fermion loops which 



wrap around the timelike direction are affected by its inclusion. The energy density e, and 
the fermion number density n, are defined by 

E = -^^^ = ^tido-foSp = -^jiTe'^M-^. - e-^M-^ (6) 
V dp V " I" ^ 2Y x,x+o x,x-0' ^ ^ ' 



n 



Here Vs is the spatial volume and V = Vs(3 the overall volume of spacetime. The final 
expression in each case is the quantity measured in the simulation, using a noisy estimator 
to calculate the matrix inverses. 

The cubic lattice has Lg lattice spacings a in spatial directions and Lt lattice spacings 
in the temporal direction. The momentum cutoff scale can be defined as A = 1/a, the 
temperature is given by T = l/{Lta) and the mass of the sigma meson is m = l/(^a), where 
^ is the correlation length of the sigma field. To reach a continuum limit one has to satisfy 
the following two conditions: A 3> T and A 3> m. The condition A ^ T requires a lattice 
with sufficiently large ^ 1. The parameters of the action should then be tuned towards 
their critical values where the correlation length ^ ^ 1. This satisfies the condition A ^ m. 



3 Monte Carlo simulations 

We performed Monte Carlo simulations of the four-fermion model in ci = 2 + 1 dimensions 
to study its phase structure in the (/x, T) plane. We used a hybrid Monte Carlo method 
described in 0, which proved to be very efficient for our purposes. Since the chiral symmetry 
is discrete we were able to simulate the model directly in the chiral limit, i.e., setting the 
fermion bare mass to zero. The three-dimensional character of the model and its relatively 
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simple form allowed us to perform simulations on large lattices and generate a large number of 
trajectories relative to state-of-the-art lattice QCD simulations. This allowed a particularly 
accurate determination of the order of the phase transition at various points in the {fi, T) 
plane. In order to check whether it is possible to detect next-to-the-leading order corrections 
to the large N calculation we set = 4. This is the smallest number of fermion species 
allowed by the hybrid Monte Carlo algorithm, because in order to have a positive semi- 
definite fermionic determinant, must be a multiple of four. 

We used the following two methods to optimize the performance of the hybrid Monte 
Carlo procedure. The first method consisted of tuning the effective number of fermion flavors 
N , which is used during the integration of the equations of motion along a microcanonical 
trajectory, so as to maximize the acceptance rate of the Monte Carlo procedure for a fixed 
microcanonical time-step dr. As the lattice size was increased, the time step dr had to be 
taken smaller and the optimal N' approached A^. For example, for the = 4 theory on a 
16 X 32^ lattice the choices dr = 0.12 and A^' = 4.026 gave acceptance rates greater than 90% 
for all couplings of interest. To maintain this acceptance rate on a 16 x 64^ lattice we used 
dr = 0.10 and A^' = 4.018. The Monte Carlo procedure was also optimized by choosing the 
trajectory length r at random from a Poisson distribution with the mean equal to f. This 
method of optimization, which guarantees ergodicity, was found to decrease autocorrelation 



times dramatically [|T3|. For most of our runs we used the average trajectory length r ~ 2.0. 
As usual, the errors were calculated by the jackknife blocking, which accounts for correlations 
in a raw data set. 



3.1 ^ = 0, r = o 

In this section we present the results of the simulations at zero temperature and zero chemical 
potential. The purpose of these simulations is to determine the scaling window of the theory 
at T = and /i = 0. The knowledge of the scaling window is necessary in order to ensure 
that the value of the lattice coupling j3 = used in the simulations at non-zero tem- 
perature and density is sufficiently close to the bulk critical coupling /3c, which means that 
lattice artifacts are negligible. In other words, we verified that another important physical 
parameter, the zero-temperature mass, mo, is sufficiently smaller than the cutoff A. We ran 
on large symmetric lattices, such as 20^ with A^ = 4 for (3 = [0.500,0.750] and determined 
the magnetic critical exponent Pmag by fitting the values of the order parameter Eq = (cr) 
to the scaling relation Sq = A{(3c - Z?)^'""". We found that (3^ ~ 0.845 and Prnag = 1.00(2) 
which is in agreement with the analytical prediction jSmag = 1 + 0(1/A^^) for the T = 0, 
= scaling According to the graph in Fig. |I|, the range of (3 for which the (So,/3) 
points fall on the 'scaling' line is /5 ~ [0.55,0.75]. This confirms that for these values of the 
coupling the lattice theory remains in the scaling window for all range of temperatures down 
to T = and effects of the lattice are negligible. 

3.2 ^^0, r = o 

In this section we discuss the results from numerical simulations of the model at non-zero 
chemical potential and zero temperature. The mean field theory approximation [^] predicts 
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that the finite /i phase transition is everywhere second order except at the isolated point 
T = 0. The issue raised by our lattice simulations in this section is whether the first order 
nature of the transition predicted by mean field theory at T = persists, when the theory 
has a small number of fermion species. In addition these simulations are motivated by our 
interest to study the difficulties in determining the order of the transition using detailed 
finite size scaling methods. We simulated the system at /5 = 0.625 because this value of P is 
deep in the broken phase and is also in the scaling window of the bulk critical point which 
means we can extract continuum physics from these simulations by using standard methods. 
The order parameter for these simulations is defined by S = (|cr|). The absolute value of a 
is necessary on a finite size lattice in order to take into account the tunneling events between 
the degenerate Z2-symmetric vacua when the parameters are tuned close to the transition. 
In the thermodynamic limit S (a). 

Figure ^ shows the order parameter S versus /z for fixed (3 = 0.625 on 16^, 22'^, 30^ and 
40^ lattices. The values of S on the bigger lattices show no fi dependence until fi reaches 
the vicinity of the critical chemical potential /ic. The variation of S becomes more abrupt 
as the volume of the system increases, indicating that a step-function type singularity will 
develop in the thermodynamic limit. The data on the 40^ lattice show a jump discontinuity 
in S as yU varies from 0.3850 to 0.3875. The fermion number density n is plotted in Fig. ^ 
as a function of /i for various lattice sizes. The transition at Hc ~ 0.3850 is seen equally well 
in this quantity. Futhermore, the variation of n near the transition becomes more abrupt 
as the volume increases, indicating again that a step-function singularity will develop in the 
thermodynamic limit. It is clear from the behavior of the various observables that relatively 
large lattices are required to simulate the theory's actual critical behavior. The discontinuous 
transition seen on a 40^ lattice is replaced by a relatively smooth curve on the 16^ lattice. 

In order to compare the fermion dynamical mass m/ at = with with the critical value 
of the chemical potential fic we calculated the fermion propagator G{'x,t) on the 40^ lattice. 
Then, we formed a zero momentum fermion propagator 

Cf{t)=Y.G{^,t). (8) 

X 

After calculating the average and the covariance matrix of Cf{t) we fitted it to the following 
functional form: 

Cf{t) = A[e™^* - (-l)*e-'^^(^*-*)]. (9) 

The value of the fermion dynamical mass is ruf = 0.403(1), which is according to our 
expectation very close to fic- The small discrepancy is due to interaction and finite size 
effects. 

Another interesting indicator which is often used to determine the order of the transition 
was proposed by Binder and Landau [jl4|, namely the reduced cumulant of the bosonic action 



Bs=l~ (10) 



S defined by 

This quantity should approach |, if the distribution of the bosonic action is described by 
a single Gaussian form with the width vanishing at infinite volume, such as in the case of 
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a second order transition; otherwise it deviates from this value. In practice we calculate 
Bs for several values of the chemical potential near the transition for each of several lattice 
sizes. Then a plot of Bs versus chemical potential is made for each lattice size. An accurate 
value of Bs,min for each lattice size is determined and then Bs,min is plotted versus y. If 
the transition is second order, this plot should have a y-intercept (representing the infinite 
volume limit) approaching the value of two-thirds. In order to locate Bs^mm accurately we 
performed simulations at very small intervals of ^ near the transition. We accumulated 
approximately 30,000 to 40,000 trajectories for the various lattice sizes at each value of 
near Hc- The cumulant as a function of is plotted in Fig. ^ for L = 16, 18, 22, 30. In 
Fig. 1^ we plot Bs^min versus inverse volume ^. The straight line in Fig. |^ represents the best 
linear fit to the data and its y-intercept is equal to 0.66475(14), which is distinctly different 
from the value of |, indicating clearly that the model has a first order finite density and 
zero temperature phase transition. A similar result was obtained from the analysis of the 
reduced cumulant behavior of the energy density. 

Our simulations were unable to distinguish a possible nuclear matter liquid-gas transition 
from the chiral transition. In the large N approximation the two transitions coincide [Q. 
The data in Fig. ^ for n versus /i show only one phase transition at /Xc- The smooth tail in 
the density before the transition is a finite size effect which dies out as the volume of the 
system increases. We extended our search for a non-zero n at /x < /Xc by performing new 
simulations on large lattices (30^, 40^, 54^) at /5 = 0.45 which is far away from the continuum 
limit. In this case the correlation length is smaller than before which implies smaller finite 
size (temperature) effects. The chiral phase transition at this value of the coupling takes 
place at /Xc = 0.692(1). Again, we did not detect a size-independent saturation in the fermion 
number density at /i < /ic. At /i = 0.690 the fermion number density n is 0.0118(4) on the 
30^ lattice, 0.0040(2) on the 40^ lattice and 0.00126(5) on the 54^ lattice. This implies that 
for = 4 the next-to-the-leading order corrections are very small if not zero, in the sense 
that either the two transitions are very close or the nuclear matter liquid-gas transition is so 
weak that is smoothed out by finite size effects or it doesn't exist in this model. 



3.3 ^^0, Tt^O 

In this section we present the results from the simulations of the four-fermion model at non- 
zero temperature and chemical potential. The purpose of these simulations is to study the 
phase structure of the model in the (/x, T) plane. The large result predicts a second order 
phase transition for T > 0. We varied T by changing either the lattice temporal extend 
Lt or the lattice spacing a which according to the discussion in Sec. |3.1| vanishes when the 
coupling [3 is tuned to its T = yU = critical value (3c ~ 0.845. 

In Fig. ^ we plot the normalized order parameter Ti{T, fi)/Ti{T = 0,fi) as a function 
of fi/rrif at different values of T. This normalization is essential in order to convert the 
order parameter and the chemical potential into physical units by getting rid of the lattice 
spacing. It therefore enables us to compare the behavior of the order parameter as a function 
of the chemical potential at different temperatures on lattices with different lattice spacings. 
According to the discussion in Sec. P]T| the lattice discretization effects cannot be neglected 



when P is smaller than 0.55. The parameters for the curves in Fig. from left to right in terms 
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of lattice size, T/TJ are: (0.65, 8 x 48^, 0.71), (0.65, 16 x 48^, 0.36), (0.55, 16 x 48^ 0.23), 
(0.55,24 X 48^,0.16), (0.625,40^,0.0). It is clear from the shapes of these curves that the 
transition becomes sharper as we decrease the temperature. In Fig. |^ we map out the phase 
diagram of the model in the {fi/mf,T/mf) plane. As the temperature increases from zero, 
the chiral condensate 'melts' at a smaller value of the critical chemical potential Hc{T) which 
suppresses energetically the quark— anti-quark pairing. 

After having determined the order of the phase transition at zero temperature (see pre- 
vious section), we extended our work to determine the order of the transition at the three 
points of Fig. ^ with the lowest non-zero values of T jT^. At these points the variation of 
E(T, /i)/S(T = 0,/i) as a function of n/nif is sharp near the transition. 

We first present the results for the order parameter E and the fermion number density 
n versus /i in Figs. |] and ^ respectively. In each case we plot the values of the observables 
for various values of the lattice spatial extend L^. For T/T^ = 0.36, = 32,48,64; for 
T/Tc = 0.23, Ls = 32,48,64 and for T/T^ = 0.16, = 48,72. These figures show a 
marked contrast with the corresponding ones at T = 0, /i 7^ discussed in the previous 
section. In each case an asymptotic behavior is approached without mutual crossings of 
curves for different Lg, in contrast to Figs. |^ and 0. However, as the temperature is lowered 
the variations of the thermodynamic observables become sharper near the transition and 
it is therefore very difficult to distinguish between a weak first order and a second order 
transition. It should also ne noted that in Fig. |^ the values of the non-zero n before the 
transition do not depend on Lg. This implies that the tail in n for fi < fic is not a finite Lg 
effect; it is instead a finite Lt or non-zero temperature effect. 

For a more accurate determination of the order of the transition we performed again a 
more detailed finite size scaling analysis. We estimated the reduced cumulants Bs for the 
bosonic action and Be for the internal energy density. For a precise calculation of these 
quantities we accumulated 30, 000 — 50, 000 trajectories for the various lattices near fic at 
T/Tc = 0.36 and 90,000 - 120,000 trajectories at T/T^ = 0.23 and T/T^ = 0.16. We 
plot Be and Bs versus fi in Figs. [10-14]. Figure |15| shows the minimum values of Be and 
Bs as functions of the inverse volume. The straight lines represent the best linear fits to 
the data. Their y-intercepts for the three different parameter sets are : Bs = 0.66667(3), 
Be = 0.66668(3) for T/T^ = 0.36; Bs = 0.66628(13), Be = 0.66641(13) for T/T^ = 0.23 
and Bs = 0.6660(3), Be = 0.6658(4) for T/Tc = 0.16. The results of this analysis confirm 
that the transition at T/Tc = 0.36 is a second order transition. They also suggest that 
at the two lower temperatures the transitions are weak first order, although the possibility 
that these might also be second order transitions cannot be excluded decisively. Studies on 
much larger lattices are required in order to detect the order of these transitions conclusively. 
Nevertheless, the steep variations of the thermodjTiamic quantities near the transitions at 
T/Tc = 0.23, 0.16 together with the deviations of the thermodynamic limit of the cumulants' 
minima from the | value support the fact that the system is close to its tri critical point. 

4 Conclusions 

In this manuscript we presented the results of Monte Carlo simulations of the three-dimensional 
four-fermion model with = 4 fermion species at non-zero temperature and density. With 
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the exception of confinement tliis model incorporates most of the essential properties of 
QCD. The great advantage of this model is that it is possible to perform lattice simulations 
at non-zero chemical potential because its fermionic determinant remains real at /z 7^ 0. We 
demonstrated using finite size scaling techniques that the /i 7^ 0, T = phase transition is a 
first order transition. Our analysis did not succeed to distinguish the nuclear matter liquid- 
gas transition from the chiral phase transition. Either the two transitions are extremely close 
or the nuclear matter liquid-gas transition which is expected to occur at yU < /i^ is very weak 
and is smoothed out by finite size effects or it doesn't exist in this model. 

We have also shown that the second order nature of the T 7^ 0, = phase transition 
remains stable for T/T^ > 0.23 and /i//ic < 0.97. This imples that the large N approximation 
(or mean field theory description of the transition) which predicts a second order phase 
transition everywhere except at the isolated point T = describes the order of the 
transition on the largest part of the critical line in the (T, /i) plane. In this sense the model 
has a "soft" behavior on the biggest section of the critical line. It should be noted that 
the mean field approximation does not predict the correct two-dimensional Ising universality 
class of the second order phase transitions line. In that case non-perturbative effects 
are crucial in the critical region and the mean field results are modified by corrections. 
However, the results of this study show that corrections to the large limit for other 
characteristics of the model, including the position of the tricritical point and the position 
of the nuclear matter liquid-gas transition, are very small indeed. 

A simple continuity argument implies that the model has a tricritical point on the critical 
line at a low temparature and very high density. Although our simulations suggest that at 
low T and high /i the phase transition is a weak first order transition, studies on much 
larger lattices would be required to "locate" the tricritical point unambiguously. Various 
approaches to QCD with two massless quarks at finite temperature and density suggest the 
existence of a tricritical point on the boundary of the phase with spontaneously broken chiral 
symmetry. The position of the tricritical point in two-flavor QCD was estimated recently 
using a random matrix model [15| and a Nambu— Jona-Lasinio model [16] as Tu- ~ 100 
MeV and fitr ~ 600 — 700 MeV. Possible experimental signatures suggested in [T^ should 
allow future experiments at RHIC and LHC to provide information about the location and 
properties of this point. 

Another conclusion drawn from our work is that lattices with very large spatial extends 
are needed in order to observe the order of the transition. Our studies on small volumes 
highlight the difficulties that must be faced when trying to understand the critical behavior 
of the thermodynamic limit by extrapolation from systems away from this limit. It was also 
shown that qualitative signatures of a first order phase transition such as the variation of 
various thermodynamic quantities near the transition can be very misleading unless they are 
accompanied by a detailed finite size scaling analysis. This should be a warning to the lattice 
community that a quantitative understanding of the critical behavior of QCD at non-zero 
chemical potential lies along a difficult road. 
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Figure 4: The bosonic action cumulant Bs vs. at T = for L = 16, 18, 22 
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Figure 5: The minimum values of the bosonic action cumulant Bs vs. inverse volume at 
T — 0. The solid line represents the best linear fit to the data. 
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Figure 6: E(T, //)/E(T — 0,fJ,) vs. ///m/ for different values of temperature T. 
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Figure 7: Phase diagram in the {T/ruf, n/mf) plane. 
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Figure 8: E vs. // for T/T^ = 0.36, (L^ = 32,48,64); T/T^ = 0.23, {Ls = 32,48,64) and 
r/Te = 0.16, (L, =48,72). 
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Figure 9: n vs. // for T/T^ = 0.36, (L, = 32,48,64); T/T^ = 0.23 (L^ = 32,48,64) and 
T/Tc = 0.16 (Ls = 48,72). The sizes of the error bars are comparable to the sizes of the 
data points. 
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Figure 10: The bosonic action cumulant Bs vs. /i for T/T^. = 0.36 and Lg = 32,40,48,64. 
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Figure 11: The energy density cumulant Be vs. for T/Tc = 0.36 and Lg = 32,40,48,64. 
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Figure 12: The bosonic action cumulant Bs vs. /i for T/T^. = 0.23 and Lg = 32,40,48,64. 
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Figure 13: The internal energy cumulant Be vs. /i for T/Tc = 0.23 and Lg = 32,40,48,64. 
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Figure 14: The bosonic action and energy density cumulants Be and Bs vs. fj, for T/Tc — 
0.16 and = 36,48,72. 
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Figure 15: The minimum values of the bosonic action and energy density cumulants vs. 
inverse volume for T/T^ = 0.16,0.23,0.36. The solid lines represent the best linear fits to 
the data. 
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